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Monte Carlo simulations using a hybrid quantum and classical mechanical potential were 
performed for crystal and amorphous- like HCl(H 2 0) n clusters (n < 24). The subsystem 
composed by HC1 and one water molecule was treated within Density Functional Theory, 
and a classical force field was used for the rest of the system. Simulations performed at 
200 K suggest that the energetic feasibility of HC1 dissociation strongly depends on its 
initial placement within the cluster. An important degree of ionization occurs only if HC1 
is incorporated into the surface. We observe that local melting does not play a crucial role 
in the ionization process. 
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1 Introduction 



Atmospheric chemistry is a research area in which many relevant processes occur in het- 
erogeneous environments, such as the surface of solid particles and within liquid droplets. 
In particular, investigations connected with the stratospheric ozone layer have proved that 
ionic solvation of HC1 at the surface of ice crystals is an important source of chlorine atoms, 
which may ultimately induce ozone-destroying chain reactions JD, ^|, [|. 

Simulations of HC1 dissociation at ice surfaces using classical force fields have recently 
been reported j|, ||. These are based on parametrizations of the potential energy surface 
which are derived from gas phase calculations for the isolated HC1(H20) dimer and the 
ionic complex Cl~ + H 3 + . Situations like those described above, in which a chemical re- 
action is strongly influenced by the environment, are rather delicate, and a purely classical 
approach risks of exhibiting problems of potentials transferability. A quantum mechanical 
semiempirical study has also been reported for HC1 solvated in water clusters ||. This 
calculation does take into account in a better way the effects of the environment, but it 
shows a very poor performance for the isolated HCl-acceptor water subsystem. This is a 
consequence of the limitations of the semiempirical description of the quantum mechani- 
cal Hamiltonian. Full ab initio Car-Parrinello simulations of HC1 dissociation have been 
performed, although in a bulk water environment fTj. 

In order to use an accurate electronic structure technique and be able to sample ade- 
quately configuration space at an affordable computational cost, we have devised a hybrid 
approach [^|, |9j in which the HCl-acceptor water subsystem is treated at the Density Func- 
tional Theory (DFT) level [TtJ and the rest of the system is modeled using the TIP4P 
potential for water We have also analyzed the role of the initial conditions and the 
local melting on the energetic feasibility of HC1 dissociation in crystal-like and amorphous- 
like clusters HC1(H 2 0)„ (n < 24) at a temperature of 200 K using Monte Carlo simulation 
techniques. 
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2 The hybrid QM/CM strategy 



The computational scheme is constructed by partitioning the system into a quantum me- 
chanical (QM) and a classical mechanical (CM) region M. Considering N c atoms in the 
classical subsystem with coordinates and partial charges {Rj, qi,i = 1,---,N C } and N q 
atoms in the QM region with coordinates and nuclear charges {r a , z a , a — 1, • ■ • , N q } the 
total energy can be written as: 

Nc ( \ Nc Nq 

E[p) = E KS [p] + / | _i | dr + EEWI Ri~ T a |) + ,J^ , ]+ E CM . (1) 



In this equation the first term is a purely quantum mechanical piece given by the standard 
Kohn-Sham expression |12|] . The electronic density p is obtained by solving a Kohn-Sham 
set of equations self-consistently, where the external potential contribution to the Kohn- 
Sham operator includes the electrostatic interaction with the CM region, as given by the 
second and third terms of expression (1). The second term accounts for the electrostatic 
interaction of the charges representing the atoms (or molecules) situated in the CM region 
with the electronic charge distribution, while the third term corresponds to the Van der 
Waals and electrostatic interactions between the nuclei in the CM region and those in the 



QM region. TIP4P parameters JXTJ] were used for O and H, and Lennard- Jones parameters 
for CI were taken from Hl3 |. The last term, Eqm, is the classical solvent contribution, 
and has been modeled with a flexible TIP4P potential for water which includes harmonic 
stretching and bending intramolecular terms extracted from extensive ab initio calcula- 
tions 0]. The electrostatic interactions between nuclei in the QM region are included in 
the Kohn-Sham expression (first term). 

For the QM region, computations are performed at the generalized gradient approximation 
(GGA) level. The correlation part is composed by the parametrization of the homogeneous 
electron gas due to Vosko [T5| and the gradient corrections given by Perdew |TjJ. The local 
exchange term was supplemented with the gradient corrections proposed by Becke [|I~7f 



The exchange-correlation contribution to the potential and the energy is calculated by a 
numerical integration scheme based on grids and quadratures also proposed by Becke |]TJ| . 
Gaussian basis sets are used for the expansion of the one-electron orbitals and also for 
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the additional auxiliary set used for expanding the electronic density. Double zeta plus 
polarization basis sets have been employed for CI, O and H |T!|. Auxiliary sets were also 
taken from Ref. 19 . 



In order to check the accuracy of the QM part of the Hamiltonian, geometry optimiza- 
tions and vibrational analysis have been performed for isolated HC1, H2O and HC1(H20). 
Structural results are shown in Table 1, together with results recently obtained at the MP2 
level ||2 Of . The agreement between DFT-GGA and MP2 results as well as with available 
experimental data is rather satisfactory. This is consistent with previous work [21], ^2[ in 
which DFT calculations at the GGA level proved to perform well for hydrogen-bonded 
dimers. 

The performance of the QM/CM approach was tested by computing the binding ener- 
gies, structural parameters and vibrational frequencies of the clusters HCl(H20) n (n=2,3), 
considering the subsystem formed by HC1 and the acceptor water molecule as the QM sub- 
system, and the remaining water molecules as the CM part. For these clusters, MP2 and 
also some experimental results are available. Selected structural parameters are shown in 
Table 1. The agreement with the MP2 computed values is again reasonable. An increase of 
the HC1 and a decrease of the OC1 bond length with cluster size can be observed, implying 
that the H-bond strength increases with the number of water molecules in the cluster. 

Binding energies and the vhci vibrational stretching frequencies for HCl(H20) n (n=2,3) 
are reported in Table 2, compared with experimental data and MP2 results. A red shift 
in the HC1 stretching frequency is experimentally observed upon complexation with water 
molecules, and reproduced by theoretical calculations. This implies that proton transfer is 
increasingly favored in larger clusters. Interaction energies for the HC1(H20) complex show 
a good agreement with MP2 calculations. Results for larger clusters show an overestimation 
of binding energies seemingly because of the use of a TIP4P classical potential parametrized 
for bulk water. However, the errors are expected to become less important for larger 
aggregates, as one approaches the bulk situation. 
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3 Monte Carlo simulations 



Finite temperature properties were simulated using a Monte Carlo (MC) technique [[TjJ . 



The MC moves consisted of random changes in the positions of all the particles simultane 
ously (including intramolecular solvent motions), with maximum displacements indepen 
dent of their respective masses. The standard Metropolis sampling algorithm was used [f28[ 



and the maximum displacements were adjusted to give an overall acceptance ratio of about 
50%. Ensemble averages were calculated over 15000 trial moves in all cases, after 4000 
moves of equilibration. All simulations were carried out at 200 K, a temperature which is 
characteristic of stratospheric conditions. 

We have considered the following situations: 



1. HC1(H20), hereafter referred as Case 1. 

2. HCl(H 2 0)i6 amorphous-like clusters. The initial conditions for these clusters have 
been obtained by running classical MC simulations at 200 K, in which the structure 
of the HC1(H20) dimer was constrained during the course of the simulation. This 
was achieved by fixing the HC1 bond length to 1.34 A(Case 2A) and 1.90 A(Case 
2B), respectively. In case 2A, the HC1 molecule remained at the periphery of the 
cluster, while in case 2B it was incorporated into the surface. 

3. HC1(H20)24 crystal-like clusters. The initial conditions were generated by isolating 
a fragment of two bilayers of hexagonal ice, composed by 25 water molecules, and 
replacing an appropriately oriented water molecule with an HC1. In the first case we 
replaced a water molecule situated in the outer layer (case 3A). In the second and 
third cases, the water molecule replaced was selected in the second layer (cases 3B 
and 3C). In case 3B, the orientation of the HC1 molecule was chosen such that it 
was H-bonded to a water acceptor molecule located in the first layer, while in case 
3C the HC1 was H-bonded to a water molecule located in the second layer. These 
are typical configurations that are likely to be found during the ice growth process 
under stratospheric conditions 0. 
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Schematic views for cases 1, 2 A, and 2B are shown in Figure 1, and for cases 3A, 3B, and 
3C in Figure 2. 

Radial distribution functions g(r) for H-Cl, H-0 (acceptor water) and Cl-0 (acceptor water) 
are presented in Figure 3 for cases 1, 2 A, and 2B and in Figure 4 for cases 3A, 3B and 3C. 
It can be observed that an important extent of ionization occurs in cases 2B, 3B and 3C, 
i.e. in those situations where the HC1 is incorporated into the surface, instead of remaining 
as an adsorbate. The degree of ionization, however, is not complete. This can be seen in 
the first peak of g(r) for H-0(acceptor water), which lies at about 1.2 - 1.3 A, while the 
optimized HO bond distance in [H 3 0] + is about 1.0 A. 

It can be observed in Figure 3 that no ionization occurs in case 2A, where the HC1 peak 
in the g(r) remains at about the equilibrium distance of the isolated HC1 molecule. The 
different behavior observed in cases 2A and 2B can be explained in terms of the solvation 
of the products, which is determined basically by the initial conditions. [H 3 0] + prefers 
trigonal coordination, and situations in which it acts as an acceptor in H-bonds are un- 
favorable. On the other hand, Cl~ prefers maximum H coordination. In case 2B, [H 3 0] + 
would be trigonally coordinated as well as the chloride ion, but in case 2A the chlorine is 
found in the periphery of the cluster and solvation is rather poor. 

Figure 4 shows that there is no dissociation in case 3A. This is because CI results with 
only coordination 2 and [H 3 0] + would act as an acceptor of an H bond (tetrahedral coor- 
dination). In both 3B and 3C cases dissociation occurs. The larger degree of ionization 
observed in Case 3B is due to the fact that, while CI always exhibits a trigonal coordina- 
tion, in case 3C the [H 3 0] + is tetrahedrally solvated, and in 3B it has the optimal trigonal 
coordination. It is also interesting to remark the different behavior observed for the g(r) 
for O-Cl in the different simulations. In the case in which HC1 is in the outer layer (case 
3A), it peaks at about 3.00 A and in cases in which it is in the second monolayer (3B 
and 3C), it peaks at 2.76 A and 2.77 A, respectively. The same trend is observed in the 
amorphous-like clusters (2 A and 2B), for which g(r) peaks at a larger value when HC1 is 
not dissociated. 
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Ensemble averages for CI Mulliken population, H-Cl, H-0 (acceptor water), Cl-0 (acceptor 
water) distances and binding energies are given in Table 3. More negative averages for 
the CI Mulliken populations are consistent with the large degree of dissociation observed 
in cases 2B, 3B and 3C. It can also be noted that larger average binding energies per 
molecule are associated with the better solvated (larger extent of ionization) situations. In 
all simulations the clusters remained solid-like, at least in the region of phase space sampled 
during our simulations. Values of Lindemann's relative rms bond length fluctuations were 
typically 0.02. Melting phenomena have not been observed, even in the simulations with 
an important degree of dissociation. 



4 Conclusions 



We conclude that the energetic feasibility for HC1 ionization in solid-like clusters strongly 
depends on the initial placement of the HC1 within the system, which in turn determines 
the solvation properties of the products. Local melting phenomena turn out not to be 
necessarily related to the dissociation process. Our results on crystal-like clusters reinforce 
the conclusions of Ref. H, in which simulations of HC1 incorporated into bulk-ice surfaces 
were performed using classical potentials. Moreover, we have shown that the same conclu- 
sion holds for amorphous-like clusters. In the case of HC1 adsorbed on top of ice surfaces it 
appears that the HC1 dissociation process would not be energetically favorable ||. These 
observations also show that the accuracy of the Hamiltonian description plays a funda- 
mental role in these studies. The QM/CM Monte Carlo scheme proposed in this work 
provides an accurate tool for modeling chemical reactions in heterogeneous environments. 

Before closing this article, we would like to make a final comment concerning ergodicity and 
proper sampling of all relevant fluctuations. During our MC runs the systems remained 
well-equilibrated and we did not observe any signature of transitions between the different 
solvation structures described in the previous paragraphs. This clearly shows the presence 
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of a high free energy barrier - in comparison to normal thermal energies - implying that, 
in principle, the feasibility of the dissociation process would be strongly dependent on the 
initial solvation conditions, i.e. on the details of the growth process. In any event, one 
would tend to believe that the more energetically favorable configuration, namely the one 
with the larger negative solvation energy (2B or 3B in our studies, see Table 3) would cor- 
respond to the most stable configuration from the thermodynamic point of view. However 
to be certain, a more complete analysis involving the computation of relative free energies 
between the different solvation structures is necessary; this would allow us to estimate not 
only equilibrium information but also information about rates of interconvertion between 
different solvation structures. Work in this direction is currently being undertaken. 
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Figure Captions 



Figure 1: 

Schematic view of initial conditions. HC1(H 2 0) ( Case 1) and HCl(H 2 0)i 5 (Cases 2A and 
2B). Only H in the QM subsystem are shown. Relevant H bonds are represented with 
dashed lines. 

Figure 2: 

Schematic view of initial conditions. HC1(H 2 0) 2 4 (Cases 3A, 3B and 3C). Only H in the 
QM subsystem are shown. Relevant H bonds are represented with dashed lines. 

Figure 3: 

H-Cl (solid line), H-0 (acceptor water) (dashed-dotted line) and Cl-0 (acceptor water) 
(dashed line) radial correlation functions, for cases 1, 2A and 2B. (distances in A) 

Figure 4: 

H-Cl (solid line), H-0 (acceptor water) (dashed-dotted line) and Cl-0 (acceptor water) 
(dashed line) radial correlation functions, for cases 3A, 3B and 3C. (distances in A) 
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TABLE 1: Selected optimized geometrical parameters for HC1, H 2 and 
HCl(H 2 0) n (n=l,3) with bond lengths in A and angles in deg. <OHCl is 
the hydrogen bond angle and dO- ■ CI and dO- ■ H the hydrogen bond lengths. 



DFT a MP2 b MP2 c Exp. 



HC1 


d HC1 




1.286 


1.271 


1.281 


1.275 d 


H 2 


d HO 




0.981 


0.961 


0.968 


0.958 d 




<HOH 




104.8 


103.5 


104.8 


104.5 d 


HC1(H 2 0) 


d HC1 




1.320 


1.287 


1.302 






d ■•• 


CI 


3.095 


3.196 


3.120 


3.2149 e 




d ■•• 


H 


1.776 


1.910 


1.818 






<0HC1 




176.6 


176.7 


178.7 




HC1(H 2 0) 2 


d HC1 




1.343 


1.303 


1.326 






d ■•• 


CI 


2.992 


3.059 


2.993 






d ■•• 


H 


1.672 


1.787 


1.688 






<0HC1 




165.6 


163.3 


166.5 




HC1(H 2 0) 3 


d HC1 




1.369 


1.323 








d ■•• 


CI 


2.923 


2.976 








d ■•• 


H 


1.558 


1.657 








<0HC1 




174.2 


174.7 







a this work. 

b 6-31g(2dp) results of Ref. 
c Poll results of Ref. 
d Ref. |23 . 
e Ref. [E4 . 
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TABLE 2: Binding energies (kJ/mol) and vhci stretching frequencies (cm 1 ) 
for HCl and HCl(H 2 0) n . (n=l,3) a 







Db 1 


Mr 2 


Mr 2 


T7 1 

Exp. 


HCl 


vhci 


2967 


3068 


2982 


2991 e 


HC1(H 2 0) 


vhci 


2512 


2841 


2709 


2659 / 












2540 9 




A E e 


23.71 


21.97 


20.57 






AE G 


14.18 


13.74 


10.93 




HC1(H 2 0) 2 


vhci 


2257 


2615 


2394 


2390 / 




A E e 


67.79 


51.27 


50.95 






AE 


50.04 


30.53 


32.86 




HC1(H 2 0) 3 


vhci 


2015 


2341 








AE e 


118.04 


89.07 








AE 


89.00 


54.70 







a AE e is the cluster dissociation energy, AE Q includes also zero point energy corrections. 
b this work. 



c 6-31g(2dp) results of Ref. |20 
d Poll results of Ref. f2| 
Ref. P 



f experimental results in Ar matrix (Ref. |25|). 
9 experimental results in N 2 matrix (Ref. PB| ) . 
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TABLE 3: Ensemble averages of CI Mulliken population, H-Cl, H-O (acceptor 
water), Cl-O (acceptor water) bond distances (A), and binding energy per 
molecule (kJ/mol). Values in parenthesis are standard deviations. 





qCl 


d H-Cl 


d H-0 


Cl-0 


E 


1 


-0.208 (0.017) 


1.314 (0.026) 


1.914 (0.118) 


3.185 (0.107) 


-7.7 (1.5) 


2A 


-0.245 (0.017) 


1.325 (0.027) 


1.752 (0.066) 


3.061 (0.057) 


-34.0 (0.6) 


2B 


-0.528 (0.036) 


1.496 (0.047) 


1.298 (0.052) 


2.789 (0.039) 


-36.6 (0.7) 


3A 


-0.280 (0.017) 


1.344 (0.030) 


1.672 (0.068) 


2.997 (0.068) 


-34.4 (0.4) 


3B 


-0.546 (0.027) 


1.510 (0.040) 


1.257 (0.037) 


2.760 (0.043) 


-35.9 (0.5) 


3C 


-0.464 (0.032) 


1.440 (0.040) 


1.335 (0.049) 


2.769 (0.043) 


-34.6 (0.6) 
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Figure 4 




